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Abstract 



We study metallization of molecular hydrogen under pressure using exact 
exchange (EXX) Kohn-Sham density-functional theory in order to avoid well- 
known underestimates of band gaps associated with standard local-density or 
generalized-gradient approximations. Compared with the standard methods, 
the EXX approach leads to considerably (1-2 eV) higher gaps and signifi- 
cant changes in the relative energies of different structures. Metallization is 
predicted to occur at a density of ~ 0.6 mol/cm^ (corresponding to a pressure 
of ~ 400 GPa), consistent with all previous measurements. 
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Despite great efforts starting with the first theoretical predictions in 1935 [|T|, the de- 
termination of the electronic and structural properties of hydrogen at high pressure is still 
extremely incomplete 0. Experimentally, it is established that hydrogen transforms to high 
pressure phases, but remains molecular up to pressures of at least ~ 200 GPa 0. MetaUiza- 
tion of solid hydrogen has been actively sought, but not yet observed, with one experimental 
team reporting no sign of metallization up to 342 GPa [|]. It is widely assumed that met- 
allization would occur either through a structural transformation to an atomic metallic 
phase, which involves dissociation of the H2 molecules, or through band overlap within the 
molecular phase itself. This latter mechanism is supported by a recent experiment-based 
equation of state ^ that, combined with Quantum Monte Carlo (QMC) calculations for 
metallic atomic hydrogen p, yields an estimate for the dissociation pressure of as much as 
620 GPa §. 

The theoretical situation is complicated by the fact that the structures at high pres- 
sures are not known, together with the well-known difficulties of quantitative predictions for 
metal-insulator transitions. Various candidate structures for the high-pressure phases (called 
"phase 11" or "BSP" below ^ 150 GPa and "phase III" or "HA" above ^ 150 GPa) have 
been proposed based on static [[f|-|T5|l, and dynamic [0,0 density-functional calculations 



and on QMC |^ investigations. Most of these structures have hexagonal and orthorhombic 
unit cells with up to four molecules. However, there are serious difficulties associated with 
the estimates of metallization pressures. The major problem is the well-known fact that the 
local-density (LDA) or generalized gradient (GGA) approximations of density-functional 
theory cause drastic underestimates of band gaps (by typically 50 - 100 %). This leads to 
much too low metallization pressures and also affects the quality of LDA and GGA total 
energies that are needed for the prediction of energetically favorable structures. Previous 



work beyond the LDA and GGA was carried out within the Xct approximation W^, the 



many-body GW framework in a first-principles and an approximate |jT5[ formulation, 
and QMC simulations 0. The first two studies were limited to a simple hep structure with 
two molecules per cell oriented along the c-axis (called "mhcp" hereafter), which has been 
found to be energetically unfavorable |l7HT0|; the GW calculations [ll5| , p!8| are not able to 
determine relative energies of structures. The QMC calculations indicated qualitatively the 
problems with the LDA calculations but did not determine gaps. 

In this Letter, we present a first-principles investigation of band-gap closure within the 
molecular phase. We employ the framework of exact exchange density-functional theory 
(EXX), which has been shown recently to yield very accurate band gaps and total energies 
for a large set of semiconductors |19[. The EXX method has crucial advantages for the 



present study. Since it treats exactly all exchange-related quantities of Kohn-Sham density- 
functional theory pO[, it is inherently self-interaction free. This remedies largely the band- 



gap underestimates that plague all LDA and GGA calculations, without an artificial band- 
gap correction and in a parameter-free way. Second, it yields band structures and total 
energies from the same calculation, which we believe to be a key requirement for the solution 
of the hydrogen problem. 

In the EXX scheme, which is explained in detail in Ref. [|l9l, total energies Etot are 
obtained from the expression 

Etot = ^0 + Eel-prot + Eh + Ex + Ec- (1) 
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Here Tq is the noninteracting kinetic energy, Eei-prot is the interaction energy between the 
electrons and the protons, Eh is the Hartree energy, E^j. the exact exchange energy, and E^. 
denotes the correlation energy, which is the only quantity that has to be approximated (the 
LDA is used in this work). Band structures {£nk} and wavefunctions (pnk for states with 
band index n and wavevector k are obtained from the Kohn-Sham equations 

(^-^ + Vprotir) + Vnir) + \4(r) + V;(r)j (/)„k(r) = e„k</>„k(r), (2) 

where V^rotlr) is the potential due to the protons, Vnir) is the Hartree potential, and 
Vc{r) = 5Ec/5n{v) with the density n(r) = 2 I0nk(r)p. The crucial part of the EXX 
scheme is the construction of the exact local exchange potential V^(r) = 6Ex/6n{r) as the 
functional derivative of the exact exchange energy with respect to the density [0. Since this 
requires repeated computation of nonlocal-exchange integrals and a linear-response function, 
an EXX calculation is much more demanding than standard LDA or GGA methods. The 
present calculations were performed using the bare proton potential for hydrogen and plane- 
wave basis sets with a kinetic energy cutoff of 36 Ry that has been employed in previous 



calculations on H2 |n],|T2|. We have also performed tests with cutoff energies of 60 Ry and 
observed that band gaps Egap and total energy differences AEtot among different structures 
were changed only by a few hundredths of an eV and a few tenths of a mRy/molecule, 
respectively. Dense k-point meshes with A^k ~ 3500/A''at special points in the Brillouin 
zone were employed {Nat denotes the number of protons in the unit cell). This guarantees 
convergence of AEtot better than 1 mRy/molecule. 

First we show how the EXX band gaps compare with LDA and GW band gaps over a 
wide range of densities, as depicted in Fig. |l| for the mhcp structure with the bond length 



and c/a ratio fixed at the values determined from LDA and extrapolations of X-ray data |1S 
(the only case for which GW calculations have been done). We can recognize several salient 
features: (i) EXX gaps are about 1.5-2 eV larger than the LDA gaps for all densities. 
Consequently, the EXX metallization density is considerably higher than for LDA. Similar 
corrections to band gaps (1-2 eV) have previously been reported for semiconductors |T9| , pT 



(ii) EXX and LDA gaps are almost linear functions of density, which holds also for the gaps 
of the other structures considered below, (iii) a linear extrapolation of the EXX data to zero 
density (isolated H2 molecules) yields a gap of 11.4 eV, close to the weighted average of the 



lowest experimental singlet and triplet excitation energies of the H2 molecule 11.6 eV 
(indicated by the left cross in Fig. |l|). 



The last point is in agreement with recent work on isolated noble-gas atoms [£3|: the 
differences between the highest occupied eigenvalue and the unoccupied EXX Kohn-Sham 
eigenvalues are very good approximations to excitation energies. This can be attributed to (i) 
the correct asymptotic — 1/r behavior of the exact exchange potential V^.(r) in Eq. (^, which 
causes the EXX spectrum of the unoccupied states to be a Rydberg series (with energies 
El < 62 < ■ ■ ■ < €00 = 0, see inset of Fig. 11]) and (ii) £00 —£homo equaling the negative of the 



ionization energy / p4|. Indeed, we find Si—Ehomo in EXX to agree very well with the lowest 
experimental excitation energies (crosses in Fig. |I|), both for the isolated hydrogen molecule 
and the molecular solid at low density. Thus, the quasiparticle band gap Egap, defined as 
the difference of the ionization energy and the electron affinity ||2^, Egap = ~ ~ Ih2-i 
differs from the lowest EXX gap by an exciton binding energy £00 — £1 @ • At high densities 
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excitonic effects are reduced, so that we expect the real quasiparticle gaps to deviate only 



slightly from the EXX gaps, just as has been demonstrated for semiconductors ||T9 . 

For densities greater than 0.35 mol/cm^, we have also carried out EXX calculations 
on other structures with hexagonal and orthorhombic unit cells (see Fig. |]) that have been 
proposed previously [0-^ as possible lowest-energy structures on the basis of LDA and GGA 
total-energy calculations. Here, the H2 molecules are tilted with respect to the z direction 
by an angle a ~ 55° and the c/a ratio is approximately 1.58 (at high pressures). In the 
structures denoted by Cmcl\, the centers of the molecules are displaced from ideal hep 
sites by a distance 8 (we normalize b such that Cmcl\ coincides with the Cmca structure). 
Proton coordinates derived from LDA calculations for these structures were used as 

input for the present EXX calculations since a complete unit-cell relaxation within the EXX 
scheme is computationally too demanding at present. 

Figure || depicts the fundamental EXX band gaps of the structures of Fig. | for den- 
sities between rii = 0.35 and n2 = 0.60 mol/cm^. [The corresponding pressures can be 
specified by our theoretical calculations or by using an extrapolated experimental equation 
of state P, p6| , p7| . We find that at the densities rii and ^2, the theoretical pressures (Pi 



and P2) are close to those of Ref. ||2^, corresponding to 100 GPa and 400 GPa; Ref. 



leads to much lower pressures at high density (Pi = 100 GPa, P2 = 325 GPa), whereas 
Ref. [27] gives higher pressures (Pi = 115 and P2 = 500 GPa).] For the Cmc2°'^, Pca2i, 
Cmc2^, and P2i/c structures, we obtain metallization densities of 0.468, 0.535, 0.537, and 
0.542 mol/cm^, respectively. Note that the use of LDA coordinates means at high pressure 
a bondlength of ro ~ 1.45 a.u. We have verified that using the experimental [r^^^^' = 1.40 
a.u.) or EXX {tq^^ ^ 1.38 a.u.) bondlength of the isolated H2 molecule [^Bj increases 
calculated band gaps by about 0.6 and 0.9 eV, respectively. For the P2i/c structure, this 
is indicated by the thin dashed lines in Fig. |^. The larger bondlength causes the predicted 
metallization density to increase up to 0.58 mol/cm^, corresponding to a calculated pressure 
of 375 GPa. 

The EXX scheme predicts that the three structures with the largest gaps (P2i/c, Cmc2\, 
and Pca2i) are the most stable ones. Our results are reported in Fig. ^ which shows the 
total energy differences among the structures for densities up to 0.62 mol/cm^. A key result 
of the utilization of the EXX functional is that the metallic Cmca structure becomes more 
stable than the insulating phases only above a density of 0.61 mol/cm^ (calculated pressure 
415 GPa). In contrast, LDA and GGA calculations [|l6l find this to be the most stable 
structure at much lower density (at quoted pressures of P ~ 140 and 180 GPa). Such a 
low metallization pressure is in severe disagreement with experiment, and we believe the 
problem is a consequence of the erroneous LDA and GGA band gaps that indirectly affect 
the total energy. The energy is decreased by populating the conduction states, an effect 
that occurs in the EXX calculations only upon much higher compression. However, the rule 
"the lower the energy, the wider the band gap" [Q is not exactly obeyed: the most stable 
structure Pca2i has only the third highest band gap. We find Pca2i to be more stable than 
P2i/c for all densities, in agreement with the LDA results of Ref. |[T^, but in disagreement 
with the LDA and GGA calculations of Refs. [|15| and [|l^] that slightly favor P2i/c. 

Zero point motion of the protons is a very difficult problem that has been the subject 
of much debate. For the present purposes there are three effects. First, the pressure is 
increased (by approximately 10% Second, the band gaps may change. The three most 
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stable structures have gaps that differ only by a few tenths of an eV. One might interpret 
this as an estimate of the influence of zero-point motion which is expected to average over 
low-energy structures. Tight-binding calculations on large cells with hydrogen molecules 
in disordered arrangements indicate effects that are similarly small. However, as the 
gaps become small, of order of the vibron energies ~ 0.5 eV, we expect the zero point 
motion to increase the gaps by a dynamic level-repulsion effect. Third, relative energies of 
different structures are changed. QMC calculations [^] and work by Straus and Ashcroft 
suggest that zero-point motion favors isotropic structures (in our case, Pca2i, P2i/c, and 
Cmc2^) with respect to anisotropic ones like Cmca. Including all these effects, we expect 
the metallization pressure to increase to ~ 400 GPa. 

Another possibility is the metallization by a structural transition to a possible monatomic 
phase. A comparison of enthalpies derived from various experimental equations of state 
with QMC calculations for hydrogen in the diamond phase yields dissociation pressures 
between 300 and 620 GPa The large range is due to the extreme sensitivity of the 

transition point to the form of the equation of state. Thus we can only conclude that 
our calculated metallization pressure is in the same range as possible transitions to other 
structures. However, this does not affect our main point that up to pressures of at least 400 
GPa, molecular hydrogen is predicted to be stable and insulating. 

In summary, we have investigated band gaps and total energies of possible candidate 
structures for compressed molecular hydrogen using a Kohn-Sham density-functional scheme 
(EXX) that treats exchange interactions exactly. EXX leads to band gaps that are 1-2 eV 
higher than in LDA (similar to gaps found recently using an approximate GW approach |I5| ) 
and, in addition, predicts changes of the relative energies of structures near the metal- 
insulator transition. In contrast to LDA and GGA calculations, the energetically preferred 
structure has Pca2i symmetry up to density 0.61 mol/cm^ (pressure ^ 400 - 450 GPa). In 
this structure there is the possibility of metallization via band overlap, which is here found to 
occur at ~ 400 GPa. Above this pressure there are three possibilities: a metallic molecular 
phase as described here; some new molecular phase that is more stable and insulating; or a 
transition to an atomic phase expected to be metallic. 

We acknowledge interesting discussions with W. Evans, A. Gorling, R. J. Hemley, J. 
Kohanoff, J. B. Krieger, P. Loubeyre, J. P. Perdew, I. F. Silvera, I. Souza, and B. Tuttle. 
This work has been supported by the Office of Naval Research under Grant No. N00014- 
98-1-0604. 
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FIGURES 



FIG. 1. Fundamental band gaps of the mhcp structure, calculated by the EXX, GW |18], and 
LDA methods. The squares represent experimental estimates for the band gap [^0|, the crosses 
denote lowest experimental excitation energies [ p^j30[| . The inset shows qualitatively the EXX 
eigenvalue spectrum in the zero-density limit. 



FIG. 2. Possible ground-state structures for solid H2 at high pressures, projected onto the xy 
plane. Pull (empty) arrows represent molecules centered on the c (c/2) plane and pointing towards 
the positive-z hemisphere. Though some of the structures have a two-molecule minimum basis, we 
have indicated the rectangular cross sections of orthorhombic four-molecule unit cells for the ease 
of comparison. 

FIG. 3. Fundamental EXX band gaps of various candidate structures for molecular hydrogen 
as a function of density. The thin long-dashed (short-dashed) line indicates the gaps of the P2i/c 
structure obtained using the EXX (experimental) bondlength of the isolated molecule. 

FIG. 4. EXX total energies of possible structures of molecular hydrogen, relative to the energy 
of the Pca2i structure. 
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